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Abstract. It has been recently pointed out that local volume fluctuations in granular packings follow remarkably well a shifted 
and rescaled Gamma distribution named the kGamma distribution [T. Aste, T. Di Matteo, Phys. Rev. E 77 (2008) 021309]. In 
this paper we confirm, extend and discuss this finding by supporting it with additional experimental and simulation data. 
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INTRODUCTION 

The description and understanding of amorphous structures is very challenging because the lack of any translational 
symmetry makes it hard to encode the structural information into a compact form. The absence of periodicity does not 
exclude however repetitiveness. Indeed, in amorphous systems several local configurations, or 'motifs', are repeated 
often. However, one must consider that these 'repetitions' typically concern similar but non identical units with small 
variations between one and the other. Furthermore, the choice of the parameters that identify the local motifs is 
somehow arbitrary and, as a consequence, depending on the detail of the description, one can gather information 
about the local structure at different levels. The identification of these motifs and the study of their variations are 
the fundamental first step toward the understanding of the structure of granular materials and of other disordered 
systems such as amorphous phases. In these disordered systems, the structure is necessarily defined in statistical terms 
and it can be characterized by the probability of occurrence of a given motif. Statistical mechanics is the theoretical 
instrument to calculate such probability. 

In this paper we focus on equal-sized bead packings both from experiments and simulations. For these systems we 
identify the local structural motifs with the Voronoi cells which are defined as the portion of space closest to a given 
bead center than to any other in the packing. In particular, we choose to consider the volume V of such cells as the only 
parameter identifying the local structural organization. The overall packing fraction <p of the sample is directly related 
to the average Voronoi volume (V) by <j) = V s / (V) with V s = nd 3 /6 the volume of a bead with diameter d. The Voronoi 
cell defines a region of pertinence around each bead. The fluctuations of the Voronoi volumes are therefore a measure 
of the local variations in the packing fraction. Indeed, the local fraction of occupied volume is <p = V s /V. In this 
paper we discuss the distribution of such volumes showing that it is very well reproduced by a kGamma distribution 
[1]. Remarkably, this functional form is retrieved in a wide set of very different systems from idealized hard sphere 
packings to glass beads in water. 

EXPERIMENTS AND NUMERICAL SIMULATIONS 
Acrylic beads in air 

The experimental data sets of bead configurations that we analyze in this paper are from the database on disordered 
packings [2] which contains structural data from experimental sphere packings obtained by X-ray Computed Tomog- 
raphy. Specifically our study concerns 6 samples (A-F) composed of acrylic beads prepared in air within a cylindrical 
container with an inner diameter of 55 mm and filled to a height of ~ 75 mm [3, 4, 5]. Samples A and C contain 
~ 150,000 beads with diameter d = 1.00 mm and polydispersity within 0.05 mm. Samples B, D-F contain <~ 35,000 



beads with diameter d = 1.59 mm and polydispersity within 0.05 mm. The two samples at lower packing fraction (A, 
B) were obtained by placing a stick in the middle of the container before pouring the beads and then slowly removing 
the stick [6]. Sample C was prepared by gently and slowly pouring the spheres into the container. Sample D was ob- 
tained by a faster pouring. In sample E, a higher packing fraction was achieved by gently tapping the container walls. 
The densest sample (F) was obtained by a combined action of gentle tapping and compression from above (with the 
upper surface left unconfined at the end of the preparation). To reduce boundary effects, the inside of the cylinder was 
roughened by randomly gluing spheres to the internal surface. The packing fraction of each of the samples is estimated 
at: A, ~ 0.586; B, ~ 0.596; C, ~ 0.619; D, - 0.626; E, - 0.630; and F, - 0.640. 



Glass beads in water 

Twelve other samples (FB12-24 and FB27) containing about 150,000 glass beads with diameters 0.25 mm are also 
analysed. The packings were prepared in water by means of a fluidised bed technique [2, 7, 8] within a vertical 
polycarbonate tube with an inner diameter of 12.8 mm and a length of 230 mm. Packing fractions between 0.56 and 
0.60 were obtained by using different flow rates with higher rates associated with lower packing fractions. After each 
flow pulse, the particles sediment forming a mechanically stable packing. 



Identification of the grain positions by X-ray computed tomography 

X-ray Computed Tomography (XCT) is used to calculate the coordinates of the bead centers. This is done by 
applying a convolution method to the three-dimensional XCT density map efficiently implemented by use of (parallel) 
Fast Fourier Transform. Furthermore a watershed method is also used to identify distinct grains. With this technique 
the estimation of the position of the centre of mass of each grain can be achieved with a precision better than 0.1% 
of the sphere diameter. This is well below the grain polydispersity that is estimated around 1 to 3 % depending on the 
sample. 



Lubachevsky-Stillinger simulations 

A set of simulated packings are produced by using a modified Lubachevsky-Stillinger (LS) algorithm [9]. The 
simulation is an event-driven Newtonian dynamics in which the spheres are considered perfectly elastic without 
any rotational degree of freedom and with no friction. The simulation is performed in a cubic box with periodic 
boundary conditions, without gravity. During the simulation, the radii of the spheres are gradually increased from a 
very loose initial state to more densely packed configurations. In these simulations the principal control parameter is 
the growth rate for the sphere radii. Small values of growth rates will result in crystallisation. To avoid crystallization 
the growth rate should be rather large, forcing the packing into "jammed" non-crystalline structures where the 
spheres cannot be further expanded at finite pressure [10, 11]. Simulations were performed by using the code at: 
http://cherrypit.princeton.edu/Packing/C-H-/ on N = 10000 spheres with initial temperature 0.1, with initial packing 
fraction 0. 1 and with a number of event per cycle equal to 20. The spheres were expanded with different growth rates 
between 2e-5 to 0.5, until a maximal reduced pressure of 10 12 was reached [12]. 



Discrete element method simulations 

We use Discrete Element Method simulation (DEM) which integrates the Newton equation of motion with both 
translational and rotational degrees of freedoms for elasto-frictional spheres under gravity [13, 14, 15]. The spheres 
interact only when overlapping, with a normal repulsive force F„ = k n ^ 2 where <^„ = d— |r,- — 7j\ is the overlap 
between grains of diameters d with centres at r ( - and rj and k„ is the stiffness parameter (k„ = d/2Y/(3(l —P 2 )), with 
Y the Young's modulus and P the Poisson ratio) [16, 17]. Tangential force under oblique loading is also considered 

as F t = —mm(\k t E,n^ 2 ^,t\, \l^F n \) ■ sign(v f ), with § ( = ff v t (t')dt' the displacement in the tangential direction that has 
taken place since the time fo when the two spheres first got in contact, where v t is the relative shear velocity and jj. 



is the kinematic friction coefficient between the spheres and k t the tangential stiffness parameter typically assumed 

2/lk n [18]. Normal visco-elastic dissipation F n = —y n ^n 2 ^ n (with 4, the normal velocity) and a viscous friction force 
F t = — YtV t [19] are also included. 

Here we report data for 64 simulations prepared by pouring, into a cylinder with a rough boundary, spherical beads 
with diameters 3mm. The cylinder had section w 22d and it was filled with 9614 beads to an height of ss 56d resulting 
in an initial packing fraction around 0.25. The beads were let sediment under gravity reaching a final mechanically 
stable state with packing fraction in a range between 0.55 to 0.64 depending mainly on the value of the friction 
coefficient (larger frictions smaller packing fractions), and also on the gravity and on the stokes constant. The time 
step has been set at 8.0e-6 sec, the grain mass is 0.003 kg and k„ =1.9e7, k, =5.6e6. The simulations were performed 
at various k s from le-3 to le2, various values of gravity from 1 to 10 and several friction coefficients from le-4 to le4. 



EQUILIBRIUM STATISTICAL MECHANICS PREDICTION FOR THE VORONOI 

VOLUME DISTRIBUTION 



Granular structures are disordered. This means that, differently form crystals, a unique "ideal" structure where all the 
grains positions are uniquely assigned does not exist. There are instead a very large number of structures that have 
equivalent global properties (packing fraction, mechanical properties, etc.) but differ in the way the grains are locally 
arranged. For these disordered packings we aim to find a relation between global functional properties and local 
structural properties and identify the probability of occurrence of specific local structural features for given global 
properties. Statistical mechanics is the theoretical framework that allows us to perform this kind of computation. 

A statistical mechanics approach for granular systems was firstly proposed by Edwards in 1989 [20]. The key 
idea is that these non-thermal systems can be described by using a formalism very similar to the one developed 
for molecular gasses by substituting the constraint on the energy with a constraint on the volume occupied by the 
system. Although this is one of the few examples of extension of classical statistical mechanics concepts to a-thermal 
systems, the Edwards' approach is rather straightforward. Any reader with some familiarity with thermal physics 
and classical statistical mechanics will recognize that the forthcoming statistical description of granular systems is 
formally identical to the one for molecular gasses with '£" substituted with 'V. However, conceptually, the approach 
is not trivial because in granular systems we lack mechanisms equivalent to temperature and molecular chaos that 
allow thermal systems to explore homogeneously the phase space. In granular systems energy is dissipated in inelastic 
collision and the system soon reaches a static state with immobile grains at mechanical equilibrium. Such state can 
only be changed by perturbing the system, injecting energy, for instance through vibrations or fluid flow [21,7]. For 
a given preparation protocol one aims to identify the probability of occurrence of some specific structural features 
and their related functional properties. In order to associate a probability to a given structure, one should (virtually) 
explore the whole set of possible structures which are attainable through a given preparation protocol and compute the 
frequency of occurrence of that specific structure within the ensemble of all attainable realizations. Within equilibrium 
statistical mechanics approach this computation is typically performed by assigning an entropy and maximizing it; 
finding in this way the configurations with maximum likelihood. It is beyond the propose and the possibility of this 
paper to fully expose the subdue issues around this kind of approach that have been debated in the literature for the 
last twenty years since 1989 [20, 22, 23, 24, 25, 26]. Here we are merely comparing the theoretical prediction from a 
statistical mechanics approach (namely Eq.4) with data from experiments and computer simulations. To have a better 
insight of our approach to this problem and for a formal deductive derivation of Eq.4, the interest reader can refer to 
[27] and references therein. Let us hereafter only briefly sketch the main ideas and the main passages of this approach. 

In analogy with the Edwards' original approach here we consider the ensemble of mechanically stable configurations 
that can be achieved by means of a given preparation protocol resulting in a given average packing fraction over a large 
number of trials [20, 22, 23, 24, 25, 26]. Here we look at the statistics of the local configurations of each Voronoi cell 
and we maximize entropy to calculate the probability p(V) for a cell with volume V in a sample with packing fraction 
where the average Voronoi volume is (V) = V s /(j). A classical equilibrium statistical mechanics theory 'a la Edwards' 
gives: 




(1) 



(2) 



v 



Here the challenge is to compute £l(V) which is the 'density of states' counting the number of microscopic configura- 
tions associated with a Voronoi volume V. 

Let us note that Eq.l is the analogous for these non-thermal systems to the Boltzmann distribution for molecular 
gasses where in this case the particle energy is substituted with the Voronoi cell volume. 



Explicit derivation of the probability distribution from a simple hypothesis 

In order to compute i2(V) here we use a very simple hypothesis: there are k 'degrees of freedom' contributing to 
the volume of each Voronoi cell. The idea is that each Voronoi cell is composed of k elements each one contributing 
independently to the cell volume V. Each of these 'elementary cells' can have arbitrary volumes larger than v m(n under 
the condition that their combination must add to a total volume V. Under this assumption £l(V) can be computed 
exactly: 

(V-kv min ) k - 1 
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(3) 



with A a constant analogous to the Debye length. Substituting into Eq.l, and by using Eq.2 we obtain for the Lagrange 
multiplier J3 = kj ((V) — kv m i„) and the probability p(V) takes the form: 

p{y) = W)m-y^f exp r(^J ' (4) 

with V min — kv min which is the minimum volume attainable for a Voronoi cell in the packing. For a packing of 
equal spheres V m ,„ is exactly known being the volume of a dodecahedral cell circumscribing the sphere which 

is V min = 5< 5 /4)/y / 2(29+13V5)rf 3 ~ 0.694J 3 [6]. Eq.4 is a Gamma distribution in the variable V - V min ; it is 
characterized by a 'shape' parameter k and a 'scale' parameter ((V) — V m i„)/k [28]. We call such a function: kGamma 
distribution [1]. Interestingly, a mathematical study for the Voronoi statistics in two dimensional point processes 
predicts a gamma distribution for the cell area distribution [29]. 
The mean of the distribution p(V) is (V) and its variance is 



g v 2 = — k '""" , (5) 
which implies 

l= M^. (6) 

This last equation gives directly the parameter k from a measure of the variance of the distribution. Therefore, there 
are no free fit parameters in Eq.4. 

It might be of some use to compute explicitly the related distribution for the local packing fraction <p which is given 
by the identity p((p)d(p = p(V)dV yielding to 

, ^_ k k (p max {(p max /(p-l)'- k ~ 1) ( h (Pmax/(p-l \ n . 

with (p max = V s /V m i n ~ 0.75 the maximum attainable local packing fraction in equal sphere packings. 



((V)-V min ) 



VORONOI VOLUME DISTRIBUTIONS FROM EXPERIMENTS AND SIMULATIONS 

We have tested the validity and resilience of the kGamma behavior over a set of several hundreds numerical simulations 
and over 18 different experiments. The simulations consisted of both Lubachevsky-Stillinger newtonian dynamics of 
frictionless hard spheres and DEM simulations of elasto-frictional spheres. The experiments include dry and fluidized 
bead samples. 

Figure 1 shows the resulting distribution of the Voronoi volumes. One can see that such distributions span a very 
broad interval of volumes with V between w l.3V s and w 2.5V S with large differences between different samples 
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FIGURE 1. (Color online), (top) Voronoi volume distributions from all the experimental and simulation data. The o refer to dry acrylic beads 
experiments and the > refer to glass beads in water. The read lines (color online) are Lubachevsky-Stillinger simulations and the green lines are 
DEM simulations. V, = 71 /6d 3 is the volume of a spherical bead, (bottom) Same plot in semi-log Y scale. 




FIGURE 2. (Color online), (top) Same Voronoi volume distributions as in Fig. 1 but plotted vs. x = (V — V,„,„)/({V) — V,„,„). (bottom) Same 
plot in semi-log Y scale. The tick full line is the kGamma function p{x) = /r(k)x k ^ i exp(— kx) for k = 13. The two dashed lines are two kGamma 
functions for k = 9 and k = 18. 
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FIGURE 3. (Color online), (left) Experimental Voronoi volumes distributions for the sample with lowest packing fraction obtained by using 
fluidized bead technique for glass beads in water (FBI 8). (right) Experimental Voronoi volumes distributions for the sample with highest packing 
fraction obtained by pouring acrylic beads in air (F). (bottom) Same plots in semi-log Y scale. The lines are the kGamma distributions with 
k = ((V) —V m i n ) 2 /Gy with G v the measured standard deviation and (V) = V s /(j> with <j> the measured packing fraction. There are no adjustable 
parameters or fits. 



shown both in the average values and in the distribution spreading. We observe that all distributions show some degree 
of asymmetry around the maximum with larger probabilities for large volume fluctuations. 

From Eq. 4, one can see that kGamma distributions characterized by similar values of k must result into similar 
behaviors when plotted as function of x = (V — V m i n )/((V) — V,,,,,,). Figure 2 shows the plot of all the distributions 
as a function of such shifted-rescaled variable. We note that all distributions collapse into a very similar functional 
form which is very well described by the kGamma function p(x) = k k /F(k)x k ^ 1 exp(-kx) with k ranging in a narrow 
interval. 

The goodness of the description of these distributions by means of the kGamma function in Eq.4 can be judged 
from Fig. 3 which shows the agreement between the prediction from the kGamma function p(V) and the measured data 
from experiments. Similar agreements are found across all samples from both experiments and simulations. It should 
be stressed that in this plot there are no adjustable parameters or fitting constants. Indeed, the only two parameters in 
Eq.4 are (V) and k which are uniquely determined respectively from the sample packing fraction ((V) = V s /(j>) and 
from the measured standard deviation of the Vorornoi volumes (k = ((V) — V m i„) 2 / Eq.6). 

The impressive fact of such an agreement is that these systems are very different (ideal Newtonian spheres, 
elasto-frictional spheres under gravity, real experimental acrylic beads in air and glass beads in water) and they are 
prepared in very different ways (pouring, tapping, fluid flows, molecular dynamics, shearing). The collapse of all 



18 



X 

$ 

X 
X 

X 
X 



X 



0.55 0.6 0.65 0/T^ 0.75 



FIGURE 4. (Color online). Behavior of k calculated from k = ({V) — V„„„) 2 /c7^ (Eq.6) vs. packing fraction <j>. The + refer to Lubachevsky- 
Stillinger simulations, the * refer to DEM simulations, the o refer to dry acrylic beads experiments and the i> refer to glass beads in water. The x 
refer instead to Lubachevsky-Stillinger simulations above (f> ~ 0.64 where a poly-crystaline phase starts to emerge. 



these distributions around a unique functional form suggests that the main driving mechanism which determines these 
fluctuations is an exchange of volume among Voronoi cells that possesses only a small number of degrees of freedom. 
Accordingly with our hypothesis, such a number is given by the parameter k. 

In Fig.4 we report the measured values of k as function of the sample packing fraction calculated from the relation 
k = ((V) —V m i n ) 2 /Oy (Eq.6). We observe a rather narrow and homogeneous range of values laying between 11 and 
15 for all samples. A sharp drop in k is observed above (j> ~ 0.64 where in the Lubachevsky-Stillinger simulations a 
crystalline phase starts to nucleate and grow. Let us note that in this phase a peak at larger packing fractions appears 
in the volume distributions (not reported in Figs.l and 2 to avoid confusion.). 

The good predictive potential of kGamma distributions for volume fluctuations in granular assemblies can also 
be inferred from the simple measure of the standard deviation. Indeed, standard deviations can be easily measured 
and require smaller statistical sets with respect to the whole distribution. Figure 5 reports the trend of the measured 
standard deviation versus the packing fraction for all the samples. It is clearly evident from the figure that they all 
follow a common decreasing trend before the crystallization onset above ~ 0.64. The line in the figure is the the 
prediction from Eq.5: a = (V s /<j> — Vmin) / vk), for k = 13. 

CONCLUSIONS 

In this paper we have shown that kGamma distributions describe very well the observed Voronoi volume distributions 
in a wide variety of systems from simulated hard spheres to real experimental packings. The agreement between the 
predicted distribution and the measures is remarkable also considering that there are no adjustable parameters or fitting 
constants. 
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FIGURE 5. (Color online). Behavior of tr„ vs. packing fraction (j>. The line is the prediction from Eq.5 (i.e. o~ = (V s /(j> — vQ)/vk) for k = 13. 
The symbols are the same as in Fig. 4. 



We have briefly explained as a statistical mechanics approach can be used to retrieve such a distribution from a very 
simple hypothesis that the Voronoi cells in the systems have k degrees of freedom associated to their volumes. Despite 
the very good agreement with the empirical data that has been shown in this paper, this hypothesis is certainly quite 
strong and unrealistic. Indeed, it might be true that in average there are k degrees of freedom per Voronoi cell but it is 
unlikely that each Vorornoi cell has exactly k degrees of freedom. A relaxation of this hypothesis on a more realistic 
ground would produce different outcomes predicting for instance other kinds of Gamma distributions. However, the 
beauty and the strength of the present approach is that we have only one parameter that is fully determined by the 
measured standard deviation. In our view, the exceptional agreement of all the experimental and numerical data with 
the prediction from kGamma distribution does not justify, at the present, the introduction of a more complicated 
theoretical framework with the consequent incorporation of extra adjustable parameters. 

Future studies will focus on the effect of polydispersity and shapes [30]. 
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